Anisotropic velocity statistics of topological defects under shear flow 
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We report numerical results on the velocity statistics of topological defects during the dynamics 
of phase ordering and non-relaxational evolution assisted by an external shear flow. We propose a 
numerically efficient tracking method for finding the position and velocity of defects, and apply it 
to vortices in a uniform field and dislocations in anisotropic stripe patterns. During relaxational 
dynamics, the distribution function of the velocity fluctuations is characterized by a dynamical 
scaling with a scaling function that has a robust algebraic tail with an inverse cube power law. This 
is characteristic to defects of codimension two, e.g. point defects in two dimensions and filaments in 
three dimensions, regardless of whether the motion is isotropic (as for vortices) or highly anisotropic 
(as for dislocations). However, the anisotropic dislocation motion leads to anisotropic statistical 
properties when the interaction between defects and their motion is influenced by the presence of 
an external shear flow transverse to the stripe orientation. 
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I. INTRODUCTION 

The small-scale dynamics of interacting defects plays 
an important role in the evolution of complex systems. In 
particular, topological defects are a common occurrence 
in systems supporting a continuous symmetry that is 
spontaneously broken in the process of a non-equilibrium 
phase transition. One central question is how the uni- 
versal properties and scaling laws near a critical phase 
transition relate to the presence and interactions of de- 
fects. The formation and evolution of topological de- 
fects is typically formulated in the framework of the 
Ginzburg-Landau theory of symmetry-breaking phase 
transitions, where defects are described as phase singu- 
larities in a complex order parameter field (rotational 
symmetry) [TJ|2]. 

The equilibrium structure of isolated topological de- 
fects is deeply rooted in their topological properties and 
is relatively well studied and understood [2]. In con- 
trast, the dynamical and statistical properties of inter- 
acting topological defects during a non-equilibrium phase 
transition are far less understood, and are the subject 
of more recent systematic analyses. Numerous studies 
have focused on the statistical properties of topological 
defect ensembles in relation to the large-scale properties 
of the system. Examples range from the quenching dy- 
namics during phase-ordering kinetics [3] 0] , the motion 
of defects in convection patterns [5], the dislocation dy- 
namics in crystal plasticity [6] [7] or the vortex filament 
motion in quantum flows [H]. A common characteristic 
of these apparently disparate systems is that they sup- 
port codimensional-two topological defects; that is dis- 
locations and vortices which in a two-dimensional space 
(2D) become point defects, and defect filaments or loops 
in a three-dimensional space (3D). One common finding 
is the presence of a robust scaling law in the local ve- 
locity statistics for these kind of defects. Recent exper- 
iments on decaying quantum turbulence in 4 He report 



that the velocity field v induced by quantized vortices 
is characterized by a v~ 3 scaling, attributed to the rare 
reconnection events between vortex filaments [9] and re- 
produced numerically in atomic Bose-Einstein conden- 
sates [TU] and counterflow turbulence [TT]. Similar ve- 
locity statistics has been observed in a discrete dislo- 
cation dynamics model of crystal plasticity [12] and in 
experiments on thermal convection in an inclined fluid 
layer [T3]. 

Theoretically, the asymptotic tail of the velocity prob- 
ability distribution P(v) can be calculated in a statis- 
tical formulation of random stationary configurations of 
point defects interacting through a logarithmic poten- 
tial in a two dimensional space (2D) [6l [14] . The model 
predicts the same tail distribution both in neutral sys- 
tems (zero net topological charge) and systems with a 
single-charge distribution. An inverse cubic scaling is 
consistent with the approximation of the nearest neigh- 
bour interaction between defects uniformly distributed in 
space [T3]. In theoretical studies of defect motion dur- 
ing phase-ordering kinetics, the inverse cubic law is re- 
lated to the annihilation events of defect loops or between 
point defects with opposite topological charges [H fTSlfTrj] . 
The coarsening during phase ordering is reflected in a 
time-dependent density of defects and their velocity dis- 
tribution F(v,t), which is characterized by a dynamical 
scaling law related to the growth law of the characteristic 
lengthscale in the ordering kinetics [13 HI] (see also [TT]). 
For non-conservative dynamics of the order parameter, 
the distribution of velocity for point defects in 2D takes 
the form F(v,t) = (v(t))~ 1 P(v/(v(t))), where the scal- 
ing function is P{x) ~ x(l + x 2 )~ 2 and the ensemble 
average velocity (v(t)) at time t is related to the average 
distance between defects Lit) at time t and scales with 
time as (v(t)) ~ 1/L(t) ~ i" 1/2 [HQS]. A different scal- 
ing exponent for the scaling function P{x) is predicted for 
defect filaments in three-dimensions (3D) [4] [15], whereas 
experiments [5] and numerics [TT] |T2] suggest the same 
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scaling as for point defects. 

In contrast to isotropic vortex dynamics, dislocation 
motion in crystals, as well as in stripe patterns, is typ- 
ically anisotropic when confined to certain gliding and 
climbing planes. In addition, dislocations often coexist 
and interact with other kinds of defects such as disin- 
clinations and grain boundaries, which makes it harder 
to study in isolation. For this reason, phase order- 
ing is much more difficult to study in isotropic stripe 
phases and polycrystalline phases, then in anisotropic 
stripes and single crystals where only dislocations are 
present (1511X5]. 

Stripe ordering is a common pattern occuring in a di- 
versity of systems from the zebra patterns to sand ripples 
and in classical fluid convection systems, where defects 
are local tears of the underlying pattern [5] . Anisotropic 
stripes or rolls develop by an uniaxial ordering of stripes 
as happens, for instance, in electrically driven convec- 
tion flows of nematic liquid crystals (electrohydrody- 
namic convection) 20J, or in thermal convection flows 
of isotropic fluids down an inclined plane [21 . 

During rclaxational dynamics, where the motion of de- 
fects is dominated by mutual interactions prior to anni- 
hilations, the statistics of the velocity components keeps 
the same form even in the presence of strongly anisotropic 
motion of defects. A numerical study of 2D phase order- 
ing after a quench from a disordered state described by a 
non- conservative time dependent real Ginzburg-Landau 
model showed that the isotropic motion of point vor- 
tices is characterized by a statistical distribution with 
an inverse cubic tail in the scaling function P(v), as 
predicted theoretically Similar statistical distribu- 

tions for the climbing (motion along the direction of the 
stripes) and gliding (motion across the stripes) velocities 
of point dislocations have been reproduced in a numeri- 
cal study of phase ordering in anisotropic stripes in two 
dimensions [23]. This is consistent with the theoretical 
understanding that the dislocation motion in anisotropic 
stripes can be in fact mapped onto a Ginzburg-Landau 
vortex-dynamics [2J. This also means that, to the lead- 
ing order approximation, the interactions between dis- 
locations are expected to be similar to those between 
vortices. 

The statistics of defect motion during non-relaxational 
evolution of the system, self-sustained or driven by an ex- 
ternal field, is less well-understood due to correlation ef- 
fects or additional driving forces apart from the mutual 
nearest neighbor interactions between defects. A self- 
sustained motion of defects is obtained in convection pat- 
terns when the mean flow due to vertical vorticity, driven 
by the undulations in the normal stripes and the presence 
of defects, acts as a self-induced drift in the motion of de- 
fects [23] . This non-trivial dynamics of defects leads to a 
spatio-temporal chaotic state also know as 'defect turbu- 
lence', which was observed experimentally in fluid con- 
vection systems [21] [25] or diffusion-reaction systems [26] . 
and studied in numerous theoretical and numerical inves- 
tigations [2~31 [2TH5U] . In this chaotic dynamical regime, 




FIG. 1. Snapshot of a measure of the charge density vector 
field for a configuration of vortex filaments in 3D simulation 
of phase ordering. The figure shows the x-component T> x of 
the Jacobian determinant defined in Sec. II. The system size 
in this simulation is 128 3 . 



a statistically stationary distribution of the number of 
defects is maintained by the defect annihilations and the 
spontaneous creations of pairs due to the phase instabil- 
ity. To leading order in the approximation of well mixed 
and independent defects, the distribution of the number 
of defects follows the Poisson statistics with mean square 
fluctuations given by the mean number of defects [27] , 
The well mixed assumption implies that defect pairs are 
being created and annihilated randomly, whereas exper- 
iments and numerical studies suggest that more often 
defects created in a pair at a given time tend to anni- 
late with each other in the same pair at a subsequent 
time |26( 130] , which means that correlations between de- 
fects are important effects in their creation/annilations 
dynamics and leads to a modified Poisson statistics in 
their number fluctuations [21] ■ A theoretical understand- 
ing for the effect of the self-induced mean flow on the col- 
lective statistical properties of defect motion is still lack- 
ing. However, experimentally measured velocity statis- 
tics during the spatio-temporal chaotic dynamics of uni- 
axial stripes in inclined layer convection are observed to 
be slightly anisotropic and the exponents in the tail dis- 
tribution of both climb and glide motion are close to 
—3 |13| . This is suggestive of a dynamical regime domi- 
nated by annihilations of dislocation pairs. 

In this paper, we consider a simpler setup where non- 
relaxational motion is driven by an externally imposed 
flow such that defects are constantly created and annihi- 
lated leading to a statistically stationary defect dynam- 
ics. This can be attained in an anisotropic stripe system 
when a shear flow is acting normal to the stripe orien- 
tation. The role of the shear flow is different from the 
commonly studied case of shear alignment of isotropic 
stripes [X3] [3XJ [32] or the buckling instability under shear 
acting along unixial stripes orientation |33j . 

The purpose of this paper is threefold: i) to present an 
efficient numerical method for tracking the position and 



velocity of topological defects and ii) its application to 
study the collective motion of dissipative vortices both 
in 2D and 3D and dislocations in 2D; hi) to report on 
numerical results where the anisotropic motion (glide and 
climb) of dislocations subjected to a simple shear flow is 
explicitly manifested in the velocity distributions, even 
though the statistics during phase ordering are similar to 
those corresponding to the isotropic motion of vortices. 

To track the position and velocity of topological de- 
fects, we implemented a numerical method inspired by 
analytical treatments of Halperin [31] and Mazenko [H] . 
The method was originally developed to locate defects 
in an 0(2)-symmetric order parameter with a Ginzburg- 
Landau relaxation dynamics in 2D. We show numeri- 
cally that this method works very well for Ginzburg- 
Landau dynamics both in 2D and 3D and it is also suit- 
able for tracking dislocations in systems controlled by 
anisotropic Swift-Hohenberg dynamics. Measuring the 
velocity statistics of vortices during relaxational dynam- 
ics, we find a universal inverse cubic tail for defects of 
the same codimension, that is point vortices in 2D and 
vortex filaments in 3D. The scaling law is directly related 
to the pairwise interactions between vortices prior to an- 
nihilation and reconnection events (in 3D). Finite size 
core effects induce a Gaussian cut-off to the v~ 3 scaling. 
A similar statistical behavior is observed in the veloc- 
ity of dislocations in anisotropic stripe patterns. Despite 
the fact that dislocations are dominated by their trans- 
verse motion, and thus are highly anisotropic, the distri- 
bution of the climb and glide velocities shows the same 
long tail behavior. In the presence of an external shear 
flow that leads to non-relaxational dynamics, the motion 
anisotropy is explicitly manifested in different statistics 
of the velocity components. While the slow motion is 
highly influenced by the shear flow, the high speed limit 
may still be dominated by the nearest neighbor interac- 
tions. 

The paper is organized as follows. Following this in- 
troduction, we proceed in Sec. II to discuss a method 
of efficiently tracking topological defects and apply it to 
a collection of vortices as well as ensembles of disloca- 
tions. Sec. Ill presents numerical results on the vortex 
velocity statistics in 2D and 3D simulations of phase or- 
dering. We discuss the statistics of dislocations during 
phase ordering and non-relaxation dynamics sustained 
by an external shear flow in Sec. IV. Concluding remarks 
and summary are provided in Sec. V. 



II. DEFECT DYNAMICS 

Here we present a numerically efficient method for 
tracking codimension-2 topological defects. The method 
is applied to Ginzburg-Landau dynamics of vortices in 
2D and 3D, as well as to Swift-Hohenberg dynamics of 
dislocations in 2D anisotropic stripes. The effect of hy- 
drodynamic interactions in the presence of an external 
shear flow is discussed in the context of dislocation dy- 
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FIG. 2. (color online) (a) Snapshot of the charge density 
field corresponding to a configuration of point vortices in 2D 
simulations. In panel (b), we show the dislocations in an un- 
derlying anisotropic stripe configuration. The lighter blobs 
correspond to defects that have a positive charge, while the 
darker blobs are the defects of the opposite charge. The sys- 
tem size in both cases is 1024 2 , while the snapshots are drawn 
from a subset. 



namics. 



A. Locating and tracking of defects 

The identification and evolution of a large population 
of defects is generally a non-trivial problem in systems de- 
scribed by continuum approaches. Moving from the field 
variables to the discrete particle variables is not straight- 
forward. In most defect studies, one resorts to various 
approximate methods to estimate the locations and ve- 
locity of defects by following their trajectories [T3"irn?ll2"2"] . 

In systems that can be described by an 0(2)-symmctric 
order parameter ip(r,t) whose evolution depicts the or- 
dering kinetics from an initially disordered state to an 
ordered state (either isotropic and homogeneous or a pe- 
riodic pattern), one can use an elegant method based on 
a transformation from the order parameter dynamics to 
the discrete defect dynamics. The analytical formulation 
of this method was pointed out first by Halperin [34] , and 
was subsequently extended by Mazenko to determine the 
velocity of various topological defects [H HH [351 |3S] • To 
our knowledge, this method has not been previously im- 
plemented numerically. We show that it is an efficient 
numerical tool used for tracking the evolution of various 
types of defects. 

The basic idea of this technique is that topological de- 
fects are located at the zeroes of the complex order pa- 
rameter field ip(r,t) |34j . The transformation from field 
to particle variables is determined by the Jacobian de- 
terminant T>(r) — \ \dil) n / drj\\, where n — 1,2 stands re- 
spectively for the real and imaginary components of the 
order parameter field, i.e. tp — ipi + iip%, and j = 1, . . . d, 
with d being the spatial dimension. Thus, for d — 2, T>(r) 
is a scalar quantity. Its sign determines the topological 
charge, i.e. q = T>(r) /\T>(r)\ = ±1, and the charge den- 
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sity is given as 

P (r, t) = 6&)V(r, t) = J2 liKr ~ n), (1) 



N 



for a collection of N point vortices. The extension to 
string defects in d — 3 is that the Jacobian determinant 
becomes a vector field T>j (r) related to the vortex filament 
density by 



Pj (r,t) = 6(i>(r,t))V j (r), 



(2) 



where the notation for the Dirac (5-function is used [3]. 

The defect velocity v is determined from the property 
of topological defects that their total charge is conserved 
(defects are created and annihilated in pairs of opposite 
charge), namely 



d t p + V • (pv) = 0, 



(3) 



with the charge density p(r,t) defined above. For exam- 
ple, in the case of point defects in 2D, the Jacobian de- 
terminant becomes V = l/(2i)CV x ip*W y ip - ^x^ v i>*), 
where tp* is the complex conjugate of ip. By differentiat- 
ing V with time, a current jW can be defined as [T5] 

^-YtW-TO), (4) 
such that the 2?-field satisfies the continuity equation 



d t V + V • J*ito = 0. 



(5) 



Summation over repeated indices is implied and e a p is 
the two dimensional antisymmetric tensor, e xx = e vy = 
and e xy = —e yx = 1- From Eqs. ^ and (|5]), the de- 
fect velocity is determined as v — /V. The defect 
velocity depends on the dynamics of the order param- 
eter through its time derivative ip(r,t). Explicitly, the 
velocity components are given by 



1 22? 
ZD 



(6) 



where %p* is the complex conjugate of -i/^-field and ip is 
the time derivative of i\) which determines the evolution 
of the order parameter. This can be generalized to d = 3, 
in which case the velocity of vortex filaments is calculated 
as [4] 



2^2 ' 



(7) 



where V 2 = Ylj=i ^j^j an< ^ the velocity vector field is 

V = (V x ,Vy,V z ). 

In the dilute defect density limit, it can be shown that 
vortex velocity defined by Eq. ([6| becomes a function of 
the phase and amplitude gradients of the order parameter 
tjj near the vortex core [35 . The formula is exact and 
applies equally well for a high density of defects. 



a= 2.92 ± 0.09i 
WOO : 
t=800 

t=1000 - 




v/<v> 



FIG. 3. (Color online) Collapsed probability distribution 
function (scaling function) of absolute velocity of point vor- 
tices in 2D during phase ordering. (Inset) PDF of the de- 
fect velocity for different core sizes (open circles correspond 
to larger core size and crosses correspond to smaller core 
size) to show that the Gaussian cut-off depends effectively 
on the vortex core size. The model parameters for the inset 
figure are A = 2.05 (open circles), A = 1.05 (crosses) and 
C = 3/20(1 + A). In the main graph, A = 1.5. Here v = |v|. 



B. Application to vortices 

Vortices are defined as the zeroes of an order parame- 
ter ip(x,t) with rotational symmetry (complex field) [2]. 
The fact that the complex field vanishes at the core of a 
defect is equivalent to a phase singularity, i.e. the phase 
of the order parameter varies discontinuously around a 
closed contour surrounding the defect. The phase 9 is 
obtained from t/j = \ip\e 10 . The shift in phase around the 
contour or the winding number, i.e. § V# • dl = 27m, de- 
fines the topological charge of the defect. A single vortex 
corresponds to a unit of topological charge, that is n = 1. 

We now consider the nonconservative evolution of a 
■0(r, t)-field described by the time dependent Ginzburg- 
Landau equation given by 



(8) 



which we simulate both in 2D and 3D. For computa- 
tional efficiency, we solve Eq. ([8| by a cell dynamical 
system (CDS) algorithm, that was originally developed 
for studying spinodal decomposition dynamics [37] and 
extensively used to study phase ordering of systems with 
continuous symmetry [3j |38l |39]. In the Appendix we 
provide a detailed description and recapitulation of the 
algorithm, for completeness, and define the parameters 
of the simulation used below. In particular, the depth 
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of the quench corresponds to the parameter A, and the 
strength of the diffusive couplings in the model are de- 
noted by C. Simulations in 2D are done on a system size 
of 1024 2 cells, while in 3D we use 128 3 cells for dx = 1. 
Unless otherwise noted, the values for the CDS parame- 
ters A, the depth of the quench, and C, the strength of 
the spatial coupling, are A = 1.5, C = 3/20(1 + A) (for 
2D), and C = 3/24(1 + A) (for 3D). Results were aver- 
aged over 48 random initial conditions, unless otherwise 
noted. 

The vortex dynamics from the Ginzburg-Landau evo- 
lution in Eq. {[sj is similar to a previous one reported 
in Ref. [25]. Here, we use a different tracking method 
for locating the defects and extend the analysis to vortex 
filaments in 3D. 

A snapshot of the charge density field for vortex fil- 
aments in 3D obtained using Halperin and Mazenko's 
method, discussed in the previous section, is shown in 
Fig. ([I]). The charge density field is directly proportional 
to the 2?-field, which is zero everywhere except along the 
vortex filaments. A similar representation is obtained for 
point vortices in 2D, where the charge field is localized 
at the vortex core and vanishes everywhere else as shown 
in Fig. ([2 (a)). Since the charge density is directly re- 
lated to the 2?-field, it means that the defect velocities 
are meaningfully defined only at the defect positions. In 
the numerical discretizations, defects are associated with 
small blobs (in 2D) or thin tubes (in 3D) with a spe- 
cific characteristic size that defines the vortex core size. 
We define the defect regions as the locations at which 
the absolute charge density is above 75% of the theo- 
retical value of | q | = 1. The values of the 2?-field are 
finite within these regions and thus the division is also 
finite. The velocity of the located defects is determined 
from Eq. Q for filaments in 3D, and Eq. ^ for point 
vortices in 2D. The time derivative ip of the ■0-field is de- 
fined by the right hand side expression in Eq. ^ , namely 

Vj = v 2 V' + V'(i-|V'l 2 )- 



C. Application to dislocations 

Here, we focus on tracking dislocations in anisotropic 
stripe patterns. A similar tracking method can be ex- 
tended to locate the defects in a crystal phase, and will 
be the subject of a separate study reported elsewhere. 

We consider the defects in a periodic pattern charac- 
terized by a preferred wavenumber k formed by stripes. 
Stripe patterns occur in a variety of systems, typical ex- 
amples being the convective rolls in Rayleigh Benard con- 
vection of isotropic fluids [5] or convective flows in the 
nematic liquid crystals [40] . in the dynamics of diblock 
copolymers [H], etc.. When the orientation of the local 
ordering is random, as in isotropic stripes, we encounter 
both isolated defects such as dislocations or disinclina- 
tions as well as grain boundaries. The coexistance of dif- 
ferent types of defects makes it difficult to analyse their 
statistics. Moreover, the ordering kinetics of isotropic 
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FIG. 4. (Color online) The scaling function of the proba- 
bility distribution function of the absolute velocity of vortex 
filaments in 3D during phase ordering. Here v = |v|. 



stripes tends to be dominated at large times by grain 
boundary slow motion, which can lead to glassy config- 
urations [18]. By fixing the orientation of the stripes 
along a preferred axis, point defects such as dislocations 
can be isolated from the other types of defects. Examples 
of anisotropic stripes are in electrohydrodynamic convec- 
tion of planary aligned liquid crystals [25], or Rayleigh- 
Bernard convection of an inclined fluid layer [13] , 

The statistics of dislocations in anisotropic stripes has 
been discussed previously by Qian and Mazenko [55] . 
where they propose a model based on an effective Swift- 
Hohenberg (SH) free energy with an additional term that 
accounts for the coupling to an external field aligning the 
stripes along a preferred direction. The stripe pattern in 
2D is represented by a real periodic field u(r,t), which 
satisfies an anisotropic SH-dynamics given by 



2^2„ 



d t u = (l-r\u\ 2 )u~ + 



(9) 



where the last term is added to impose a preferred orien- 
tation of the stripes along the vertical y-axis with c > 
being the coupling strength to the external field. The 
quench depth r > is interpreted in the context of con- 
vection patterns as the deviation from the onset of con- 
vection, r f=a R/R c — 1, where R is the Rayleigh number 
and R c is the critical R at the onset [5]. The anisotropic 
preferred orientation can be seen by linearizing Eq. ^ 
around the mode solution u ~ exp(wt+ik x x-\-ik y y) with 
the growth rate obtained from Eq. ^ as 



u = {l-r)-(\-k 2 x -kl) 2 +ckl, 



(10) 



and imposing the condition Lu(r;k x ,k y ) = at the onset 
of instability with respect to a mode of wavenumbers k x 
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and k y . The condition is found by minimizing uj with 
respect to k x and ky, i.e. duj/dk x = and duj/dk y = 0. 
This leads to k y = and k x = \/l + c/2 for c > 0. 

We consider the amplitude formulation of Eq. (J9j) with 
an additional contribution due to an external shear flow. 
We impose a shear flow that is normal to the main orien- 
tation of the stripes to allow for the nucleation of defects 
due to wavenumber shifts by shear deformation. In or- 
der to track dislocations in a complex order parameter 
field, we write the periodic field in terms of its complex 
envelope field ip(r,t), namely u(r) = y / r0(r)e lfc ' r + c.c 
Without loss of generality, we consider stripes with the 
wavevector parallel to the horizontal x-axis, i.e. k = 
(ko,0). The complex ■(/'-field satisfies an amplitude equa- 
tion derived from Eq. ^ and given to the leading order 
in r as 

a t v + %CM = Ki - M 2 )^ - £ 2 m - c4M,(n) 

where C = (V 2 + 2ik ■ V) is derived from (I + V 2 ) and 
C x = Va, + ikg comes from the gradient V x . The advec- 
tion term is determined by a velocity field, which hereby 
is taken as a simple shear flow 7 = voyx, and the last 
term is added to impose a preferred orientation of the 
stripes long the vertical y-axis. 

We integrate numerically Eq. (ll| ) using a 4 t?l -order 
Runge-Kutta scheme and a spherical approximation for 
the gradients |42| (see Appendix) on a square domain 
of size 1024cfe x 1024rfx. The time step is dt = 0.05 
and the spatial resolution is dx — 7r/4 so that about 8 
grid points are used to resolve the pattern wavelength 
A = 27r/fc , with k = 1. The other parameters are set 
to c = 1, r — 1 and v$ is a changing parameter. In the 
absence of shear, period boundary conditions on all sides 
are used. At a finite shear rate, we impose a zero flux 
boundary conditions of the upper and lower boundaries 
and periodic conditions on the lateral boundaries. 

Dislocations are efficiently located as the zeros of the 
complex envelope field ip using Mazenko's algorithm. 
In Fig. |(2] panel (b)), we illustrate a stripe configura- 
tion with the location of dislocations and their topo- 
logical charge proportional to the Jacobian determinant 
||cfy/? n /0r/||. The velocity of dislocations is obtained 
using Eq. ([6| with the evolution of the order param- 
eter given by the right hand side of Eq. (11), namely 
i> = r(l-W\ 2 )i>-C 2 ty]-c£ 2 M. 



III. VORTEX STATISTICS 

To determine the velocity statistics of vortices, we ini- 
tiate the system in a disordered state and follow the or- 
dering kinetics dominated by the initial formation and 
subsequent coarsening of topological defects. At a partic- 
ular time, we calculate the local defect velocities v using 
Mazenko's method as described above. We save the ab- 
solute values, v = \v\, every few time iterations and run 
the system from 48 random initial conditions. This way, 
we compute the probability distribution function of the 
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FIG. 5. (Color online) a) Density of defects pd as a function 
of time t for different values of the imposed shear velocity 
vq. b) Data collapse of the rescaled density with the mean 



density in the statistical steady state po 
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where to 



is the cross-over time to the steady state. Inset figure shows 
how to scales with vo as to ~ 1/vo- 



defect velocity at a given time, i.e. F(v, t). In the asymp- 
totic limit t — v 00 of the coarsening dynamics, we expect 
scale invariance of the typical coarsening lenghscale, i.e. 
L(t) ~ i 1 / 2 for non-conservative dynamics (apart from 
logarithmic corrections in 2D). Hence the typical veloc- 
ity obtained as 1/L(t) scales with time as (v(t)) ~ t -1 / 2 
and corresponds to the velocity of defects in a pair prior 
to annihilation and separated by a distance of the order 
of L(t). In the simulations, (v(t)) is the ensemble aver- 
age velocity at a particular time, and when calculated 
over long times it converges to the expected asymptotic 
scaling. This dynamical scaling of the mean velocity im- 
plies also a scaling with time of F(v,t). We notice that 
the time dependence in the PDF's can be eliminated by 
rescaling the velocity variables by their ensemble average 
values at a given time, i.e. v = v/(v(t)). Analytically, 
this corresponds to the rescaling F(v,t) — t 1 / 2 P{vt 1 / 2 ). 

The scaling function P{x) of the velocity distribution 
is a function of the rescaled velocity field v/ (v(t)) and has 
a broad tail with an inverse cubic decay. This is shown 
in Fig. (J3J) for 2D simulations and Fig. Q for 3D dy- 
namics. The v~ 3 tail corresponds to the regime of large 
velocities obtained in pair interactions prior to annihila- 
tion events or, for 3D, also reconnections events. The —3 
scaling exponent is determined by the logarithmic mu- 
tual interaction potential as shown in e.g. Refs. [HI [43]. 
Since point defects in 2D and filaments in 3D are both 
codimensional-2 topological defects with the same type of 
interactions, we expect a similar scaling behavior. This 
is consistent with other numerical studies that also show 
that the tail of P(v) distribution is dominated by the v -3 
scaling both in 2D and 3D simulations [TTJ [T^] . We pro- 
vide additional numerical evidence for this scaling regime 
also during phase ordering kinetics. 

From Figs. ^ and Q, we notice that the vortex core 
structure has a drastic influence on the defect velocity 
statistics [43]. This effect is given by a Gaussian tail 
which takes over the v~ 3 regimes at very large fluctua- 
tions. The vortex core size is more pronounced in 3D 
simulations, due to numerical bounds on higher spatial 
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FIG. 6. (Color online) a) Evolution with time of the ensemble 
average of the velocity component v x = (\v x (t)\) for different 
values of the shear rate, b) Data collapse of the rescaled v x 
as a function of the rescaled time t/to. 



FIG. 7. (Color online) a) Evolution with time of the ensemble 
average of the velocity component v y = (\v y (t)\) for different 
values of the shear rate, b) Rescaled v y plotted as a function 
of the rescaled time t/to. 



resolutions. It is easier to vary the aspect ratio between 
system size and vortex core size in 2D simulations to ob- 
serve the vortex core effect. In CDS simulations, the core 
size is fixed by the parameter A. We consider two dif- 
ferent values of A, namely A = 1.05 for small cores and 
A = 2.05 for larger cores, and averaged over 2000 random 
initial conditions. The system size was reduced to 128 2 
cells. The dependence of the Gaussian-cutoff on the core 
size is shown in the inset of Fig. ([3]), while the velocity 
distribution for an intermediate core size, parametrized 
by A = 1.5, is plotted in the main graph of Fig. using 
1024 2 cells and averaging over 48 initial conditions. 



IV. DISLOCATION STATISTICS 

Next we discuss the statistics of dislocations in 
anisotropic stripes during phase ordering and non- 
relaxational dynamics assisted by a shear flow. 



A. Dynamical scaling regimes 

In the early stages, the evolution of defects is dom- 
inated by pairwise interactions leading to annihilations 
and a decrease in the density of defects. We measure the 
density of defects Pd{t) as the ratio between the effec- 
tive area occupied by the defects and the total system 
area. Alternatively, the number of defects N can be es- 
timated as the area occupied by the defects divided by 
the approximate core area of a defect. In simulations, we 
choose a large system size of 1024 2 and calculate Pd(t) 
as a function of time, starting from a random initial con- 
dition. During the relaxation dynamics, the defect den- 
sity for a very large system is equivalent to the averaged 
density over many initial conditions for a smaller system 
size, and is a smooth function of time [U] . As expected 
in the coarsening regime, the density obeys a power law 
in time as Pd{t) ~ \og(t)/t like the density of vortices 
in Ginzburg-Landau theory [22] . This behavior is shown 
in Fig. (|5^) by the data in open circles and is consis- 



tent with formal arguments that the anisotropic Swift- 
Hohenberg dynamics can be mapped onto the isotropic 
Ginzburg-Landau dynamics j2]. In Fig. ((5^,), we also 
plot the density Pd{t) as a function of time for various 
values of the shear rate vq. We notice that in the late 
stages, when long-range hydrodynamic interactions set 
in, the density of defects ceases to decrease monotoni- 
cally, and approaches instead a statistically steady state. 
In this steady state , the defect density fluctuates in time 
about a mean value because of the sporadic pair creations 
and subsequent annihilations of defect pairs. Averages 
over initial conditions would lead to a constant mean 
density po in this steady state, i.e. (pd(t))ic —> Po as 
t — > co. We also observe that the mean number of de- 
fects (N) in the steady state appears to increase with 
the applied shear v , suggesting that the creation rate of 
defect pairs depends on vq. Also, the mean square fluc- 
tuations (A 2 ) — (A) 2 in the number of defects increases 
monotonically with (N) with significant deviations from 
what would be expected in a Poisson process. It appears 
that the number statistics of defects behaves in a simi- 
lar manner to that in defect turbulence [50] ; although a 
detailed analysis of this suggestion would be beyond the 
scope of this paper. 

From the data presented in Fig. ^p), we can deter- 
mine the cross-over time £ to the statistical steady state, 
and we find that it increases to a first approximation as 
to ~ Vq 1 , as shown in the inset of Fig. ^p). It turns 
out that phase ordering assisted by hydrodynamic ef- 
fects slows down the growth of the typical distance be- 
tween defects with corrections that follow a L(t) ~ t 1 / 3 
law until saturating to the steady state. This implies 
that the steady state defect density can be estimated as 
po ~ L(io)~ 2 ~ t Q 2 ^ 3 . Equivalently, the mean density 
in the steady state increases with the applied shear as 
Po ~ v q^ 3 - As a nrs t step to see whether there is any 
data collapse associated with this scaling behavior, we 
plot the rescaled p/po versus the rescaled t/to as shown 
in Fig. <|5Jd). 

Because the normal stripes are aligned along the y- 
direction, the climb and glide motions of the dislocations 



8 



correspond respectively to the vertical and horizontal ve- 
locities. The motion is anisotropic and dominated by the 
transverse (gliding) dynamics of dislocations. We ob- 
serve that climb motion typically occurs when two dis- 
locations of opposite topological charge approach each 
other to annihilate, otherwise dislocations would move 
by gliding across the stripes. Nevertheless, when the 
two motions are driven by the local pairwise interac- 
tions between dislocations, they exhibit similar charac- 
teristics. In the early states, where the dynamics is con- 
trolled mainly by the phase ordering and annihilations 
of dislocations, the ensemble average of the net veloci- 
ties (absolute values), (\v x \}(t) and (\v y \)(t), scale with 
time as i" 1 / 2 until hydrodynamic effects set in and the 
defects are accelerating until they cross over at t to a 
statistically stationary state. This behavior is shown in 
Fig. (|6h) for (\v x \) and in Fig. ([7^) for (\v y \) correspond- 
ing toaifferent values of the imposed shear velocity vo. In 
the later stages, hydrodynamic effects become important 
and there is a transient regime where the mean velocities 
increase almost linearly with time until to after which a 
statistically stationary state is reached. Since, we run 
large scale simulations for a given realization without av- 
eraging over initial conditions, the ensemble-averaged ve- 
locities, and (\v y \)(t), correspond to time series 
in the statistical steady state. Averaging these fluctua- 
tions over many initial conditions would smooth out the 
late-stage time dependence to a constant value. 

The mean square fluctuations in the steady state are 
an increasing function of the applied shear. We observe 
that the mean value of climb velocity is approximately 
an order of magnitude smaller than the typical velocity 
of gliding. A rescaling of velocity components as a func- 
tion of the rescaled time in the units of to is presented in 
Figs. (JBJd) and ([7]b) , which however gives a poor data col- 
lapse. Improving it turned out to be a challenging task, 
one of the reasons being that there is an additional char- 
acteristic timescale given by the cross-over from the re- 
laxation dynamics to the transient period of acceleration 
prior to the steady state. It may be that this timescale 
also plays a role in the scaling function, but we have not 
succeeded in including it to our satisfaction. This is an 
unresolved issue that deserves a separate detailed study. 

The statistically stationary regime is characterized by 
fluctuations in the density of defects and their velocities 
around a mean value that depends on the imposed shear 
flow. Fluctuations in the defect density are attributed 
to sudden nucleation of dislocation pairs and their sub- 
sequent annihilation, either with the same pair member, 
or with dislocations from another pair. 

Physically, the nucleation mechanism is related to the 
phase shifts induced by the transverse shear deforma- 
tions. The reason for this is that the action of an exter- 
nal shear flow is cumulative in the phase 9 of the complex 
envelope field ip = |^|e lS and its gradients V0, similar to 
the effect of the self- induced mean flow [45]. The shear 
flow advects the underlying stripe pattern together with 
its defects. This motion induces small undulations along 




FIG. 8. (Color online) The time independent scaling func- 
tion of the probability distribution corresponding to climb 
(crosses), respectively glide (open circles) velocities during the 
phase ordering in the absense of external field. The rescaled 
velocity variable is given as Vj = Vj/{\vj\) for j = x, y. In the 
inset figure it is shown the log-log plot of the PDF calculated 
with logarithmic binning. 



the stripes which build up stresses and create distortions 
in the pattern. These distortions, which are localized into 
transverse 'shear bands' , grow with time up to the point 
where they locally tear apart the stripes releasing pairs of 
defects. The shear flow affects both the isolated motion 
of defects and, more importantly, the interaction between 
defects. The effect of the large-scale flow on the defect 
interactions becomes stronger where the defect motion is 
slower [45] . 



B. Velocity statistics 

In the absence of shear flow, the motion of dislocations 
is symmetric. Although the mean velocity in absolute 
value is non-zero and related to the mutual interaction 
forces, the velocity of dislocations averages out to zero. 
This is equivalent to the symmetric probability distribu- 
tions of the velocity fluctuations as shown in Fig. (JsJ) . 
The actual distribution is time dependent due to quench 
dynamics by annihilations. However, using the dynam- 
ical scaling behavior of the probability distribution we 
can remove the time dependence by effectively rescaling 
the dislocation velocity by the ensemble average of the 
absolute velocity at a given time, i.e. Vi = Vi / (\i)i(t)\) , 
with i = x,y for the two components. The distribution of 
these rescaled velocities corresponds to the scaling func- 
tion of the time dependent probability distribution as 
discussed previously in the context of vortex dynamics. 



From Fig. (J8|, we notice that the probability distribu- 
tions of the climb and glide velocities retain a similar 
form that is characterized by a long tail with a —3 power 
law as in the relaxational dynamics of vortices. This is 
consistent with previous numerical studies of the velocity 
statistics of defects in the anisotropic Swift-Hohenberg 
dynamics from Ref. [23] . 

At a finite shear rate and in the late stage of statis- 
tically stationary regime, the motion of dislocations is 
influenced by the imposed flow. Dislocations of opposite 
topological charges tend to move in opposite directions, 
with an asymmetry in their mean transverse motion that 
is related to the shear rate, i.e. 7 ~ ((t>+) - {v~}), while 
the climb motion remains almost symmetric. The ve- 
locity probability distribution for positive, respectively 
negative dislocations becomes the same when we rescale 
their corresponding absolute velocities as v s = |w s |/(|f s |), 
where s = ±, so that P + (v + ) = P~(v~). 

In Fig. ([9]), we plot the probability distribution func- 
tions of the rescaled gliding velocities and climbing veloc- 
ities in the steady state regime. For the transverse mo- 
tion, the small scale velocity fluctuations are normally 
distributed around the mean flow. Large fluctuations 
against the flow are due to events where pairs of dislo- 
cations of opposite charge, gliding opposite to their drift 
flow, are attracting and annihilating. These events con- 
tribute to the long left tail for P + (v x ) and right tail for 
P~(v x ). The statistics of small climb velocities are also 
influenced by the imposed shear flow and given by the 
slight asymmetry in the P s {v y ). However, the large fluc- 
tuations in the longitudinal motion are due to pair in- 
teractions prior to annihilations. These fluctuations are 
captured by the long tails in the P{y v ) distribution and 
they seem to follow the inverse cubic law, but with less 
accuracy than in the relaxational case. Typically, a nu- 
cleation event leads to a burst in the local density of dis- 
locations which will annihilate subsequently by the fast 
climbing motions. However, these large velocity events 
occur on a longer timescale than during relaxational dy- 
namics, so that the system needs to be followed longer 
in the steady state. This is computationally challenging, 
because of the prior transient acceleration period whose 
length increases as the driving force is decreased. 



V. CONCLUSIONS: 

In summary, our numerical simulations suggest that 
the velocity statistics of codimensional-two defects ex- 
hibits a dynamical scale invariance with a scaling func- 
tion that has a universal inverse cubic tail when the de- 
fect dynamics is driven by mutual pair interactions lead- 
ing to annihilations. This is valid both for point defects 
in 2D and defect filaments in 3D during phase ordering 
kinetics. Finite size effects of the defect core introduce 
a Gaussian cut-off to the v~ 3 scaling. A similar statis- 
tical behavior is observed in the velocity of dislocations 
in anisotropic stripe patterns during phase ordering. Al- 
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FIG. 9. (Color online) (a) PDF of the gliding velocity for 
positive (open diamonds and circles) and negative (open tri- 
angles and crosses) dislocations in the statistically stationary 
state at a different shear rates, (b) Probability distribution 
function of the climb velocity for positive (open circles and 
diamonds) and negative (open crosses and triangles) disloca- 
tions in the statistically stationary state. In both cases, the 
rescaled velocity variable is Vj = Vj/{\vj\) for j = x,y. Inset 
figure shows the log-log plot of the tail PDF with logarithmic 
binning. 



though the motion is highly anisotropic and dominated 
by the gliding of dislocations, the distributions of the 
climb and glide velocity fluctuations exhibit the same al- 
gebraic tail when the motion is driven by the local in- 
teractions between dislocations. During non-relaxational 
dynamics assisted by an external transverse shear flow, 
small velocity fluctuations are influenced by the mean 
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flow, whereas the asymptotically large fluctuations are 
still due to pairwise interactions. In statistically station- 
ary dynamics, anisotropic defect motion is manifested 
also in anisotropic statistics of the glide and climb ve- 
locities. In this study, we have neglected the effect of a 
self-induced mean flow in the defect dynamics. This ef- 
fect is however important to capture the spatio-temporal 
chaotic dynamics as seen from experiments. It is thus 
interesting to study in future investigations the statisti- 
cal properties of interacting defects when the combined 
effect of a large scale flow and a self-induced mean flow 
is taken into account. 
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Appendix: Calculation of derivatives in cell 
dynamical systems models 



We solve the time dependent real Ginzburg-Landau 
dynamics from E q. j8 | using the cell dynamical systems 
approach (CDS) |38 ] l39 | [46] . Numerical efficiency be- 
comes particularly important in 3D simulations and this 
method is tailored to that. The complex variable 4>(r,t) 



is replaced by a ip^j (or i n 3D) defined on a square 
lattice of size N x N (a cube of size N x N x N in 3D) 
at time n. The idea of the CDS method is to construct 
a discrete set of maps for each lattice cell such that the 
flow properties of the continuous dynamical system are 
preserved. A cell dynamics is defined by two steps: a 
local update 



4> 



(n+l) 



^/l+^(n)(A 2 -1)' 



(A.l) 



where A > 1 is a parameter that determines the global 
rate of convergence to the fixed points of the local double- 
well potential, and a global update taking into account 



the interactions between neighbouring cells 



(A.2) 



where C is a constant proportional to the phenomenolog- 
ical diffusion constant. The isotropy of the order param- 
eter being simulated naturally mandates the isotropy of 
the difference operators used to implement the coupled 
maps (see Tomita [12] )■ Oono and Puri [15] chose a 9 
point stencil to implement a "Laplacian" operator, which 
is highly isotropic for a 2D square lattice [17] . This sten- 
cil reads 



3 

dx 2 



\ NN 



1 

12 



- ^ j , (A.3) 



NNN 



where NN stands for the nearest neighbours in the dis- 
cretized lattice and NNN are the next-to-nearest neigh- 
bours for each node in the lattice. 

Considering the same isotropy requirements, the dis- 
cretization of the 3D Laplace operator reads as [39] 



V 2 / 



3 
dx 2 



&5> 

\ NN 



(A.4) 



NNN 



Since the calculation of the position and velocities of 
defects involves first order spatial derivatives of the ip 
field, an isotropic discretization of the gradients is impor- 
tant in order to reduce the underlying lattice anisotropic 
effects. We use the isotropic version of the gradients both 
in 2D and 3D. Following the idea that these operators 
have to be accurately represented in Fourier space |47j . 
we use the following stencil for 2D 



■2^ 
2ip itj -i - ijji 



- V'i-ij+i 
•ij-i) > 



(A.5) 



where i and j are the lattice indices for the x and y 
directions, respectively. Swapping indices we can obtain 
the corresponding expression for ^/ y ip. We note that this 
expression looks very similar to a 4 point first derivative 
in a 2D square lattice [IB]. For the gradients in 3D, the 
stencil reads 



V x -0 = ^-j^ Wi+l,j+l,k — 1pi-l,j+l,k + ^i+l,j-l,fc — i>i-l,j-l,k + ^i+l,j,k+l ~ 1pi-l,j,k+l + 1pi+l,j,k-l ~ 1pi-lJ,k-l) + 



+ 73- (tpi,j+l,k - i>i,j-l,k + *Pi,j,k+l - "4>i,j,k-l) ■ 



(A.6) 



with the corresponding index swap to obtain V y ip and V 2 '0. 
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